https://felixml.github.io/ECS530_h20/
Salmon Farm in Northern Norway
A large farm may contain as many fish as the entire population of wild salmon in Norway.
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(tmap)
M <- tm_shape(Nor) + tm_polygons(, col= "grey")
P <- tm_shape(Akvaf) + tm_dots(col= "red" , size = 0.001 )
P2 <- tm_shape(Akvaf) + tm_dots(col= "red" , size = 0.01 )
M + P
With about 1,000 farms, the industry produces X1000 more than the natural carrying capacity of the ecosystem.
The biomass density of fish in the sea enables a parasite to thrive. The sea lice attaches to the skin of salmon and eats it, causing damages and infection. For the industry this is very costly. It also affects wild smolt migrating down the rivers into the sea.
The government has introduce several legislation which limits the growth of the industry.
library(raster)
st_bbox(Nor)
## xmin ymin xmax ymax
## 3.904688 57.959026 31.161980 71.181252
rasta <- raster(Nor, st_bbox(Nor))
R1 <- rasterize(Akvaf,rasta, field = "LOK_NR", fun="count")
qtm(R1) + M
library(tmaptools)
Bomlo <- bb("Bomlo")
Bomlo2 <- tm_shape(Akvaf, bbox=Bomlo) + tm_dots(col= "red" , size = 0.1 ) +
M + tm_layout(bg.color="blue") +
tm_compass(type="8star", position = "left", color.dark = "white", color.light = "black", text.color="red" )
Bomlo2
library(magrittr)
##
## Attaching package: 'magrittr'
## The following object is masked from 'package:raster':
##
## extract
library(leaflet)
#Map Hardanger
tmap_leaflet(P2) %>% setView(lat=60.0, lng=5.8, zoom= 8)
#Map Bomlø
tmap_leaflet(P2) %>% setView(lat=59.8, lng=5.5, zoom= 9)
#Nearest Neighbor 1 : spdep
suppressPackageStartupMessages(library(spdep))
dnb <- dnearneigh(Akvaf, 0, 25, row.names=row.names(Akvaf))
dists <- nbdists(dnb, st_coordinates(Akvaf))
ozpo <- set.ZeroPolicyOption(TRUE)
lw <- nb2listw(dnb, glist=dists, style="B")
## Warning in nb2listw(dnb, glist = dists, style = "B"): zero sum general weights
#Nearest Neighbor 1 : nngeo
library(nngeo)
nn = st_nn(Akvaf, Akvaf, k = 3, progress = FALSE)
## lon-lat points
st_join(Akvaf, Akvaf, join = st_nn, k = 3, progress = FALSE)
## lon-lat points
## Simple feature collection with 3171 features and 2 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 4.660783 ymin: 58.0243 xmax: 30.41533 ymax: 71.01778
## geographic CRS: WGS 84
## First 10 features:
## LOK_NR.x LOK_NR.y geometry
## 157 12832 12832 POINT (8.26705 58.14398)
## 157.1 12832 14016 POINT (8.26705 58.14398)
## 157.2 12832 29957 POINT (8.26705 58.14398)
## 159 14016 14016 POINT (8.2553 58.14027)
## 159.1 14016 12832 POINT (8.2553 58.14027)
## 159.2 14016 29957 POINT (8.2553 58.14027)
## 163 29957 29957 POINT (7.06875 58.05583)
## 163.1 29957 30657 POINT (7.06875 58.05583)
## 163.2 29957 33577 POINT (7.06875 58.05583)
## 164 33577 33577 POINT (6.885683 58.06008)
l = st_connect(Akvaf, Akvaf, ids = nn, progress = FALSE)
plot(l, col = NA) # For setting the extent
plot(st_geometry(Akvaf), col = "darkgrey", add = TRUE)
#plot(st_geometry(Akvaf), col = "red", add = TRUE)
plot(l, add = TRUE)
By Sea Distance
#Buffers( What are the units here?)
buff25 <- st_buffer(Akvaf, 25)
## Warning in st_buffer.sfc(st_geometry(x), dist, nQuadSegs, endCapStyle =
## endCapStyle, : st_buffer does not correctly buffer longitude/latitude data
## dist is assumed to be in decimal degrees (arc_degrees).
plot(buff25)
## Warning in st_is_longlat(x): bounding box has potentially an invalid value range
## for longlat data
## Warning in st_is_longlat(x): bounding box has potentially an invalid value range
## for longlat data
https://www.barentswatch.no/nedlasting/fishhealth/lice
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
##
## intersect, select, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
summary(lus30$Lice)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.0000 0.0000 0.0600 0.1208 0.1900 2.1200 488
nrow(lus30)
## [1] 1064
tapply(lus30$Lice, lus30$Fylke, summary)
## $Agder
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.00000 0.03000 0.06000 0.06333 0.09750 0.13000 7
##
## $`Møre og Romsdal`
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.00000 0.00000 0.02000 0.06188 0.09000 0.36000 36
##
## $Nordland
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.00000 0.00000 0.02000 0.09016 0.13000 0.78000 104
##
## $Rogaland
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.0000 0.0225 0.1650 0.1947 0.3375 0.5000 35
##
## $`Trøndelag`
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.0000 0.0100 0.1000 0.1771 0.2500 2.1200 74
##
## $`Troms og Finnmark`
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.0000 0.0000 0.0200 0.0522 0.0500 0.4900 102
##
## $Vestland
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.0000 0.0325 0.1300 0.1612 0.2700 0.4900 130
#lm(lus30$Lice ~ lus30$Temp + I(lus30$Temp^2))
vest <- lus30[lus30$Fylke=="Vestland",]
nrow(vest)
## [1] 292
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 162 0.16 0.14 0.13 0.15 0.15 0 0.49 0.49 0.62 -0.78 0.01
library(colorspace)
##
## Attaching package: 'colorspace'
## The following object is masked from 'package:raster':
##
## RGB
pal <- heat_hcl(10)
mvest <- st_as_sf(vest, coords = c("Lon", "Lat"), crs = 4326)
plot(mvest["Lice"], nbreaks=7, breaks="fisher", col=pal)
## Warning in plot.sf(mvest["Lice"], nbreaks = 7, breaks = "fisher", col = pal):
## col is not of length 1 or nrow(x): colors will be recycled; use pal to specify a
## color palette
library(spatstat)
## Loading required package: spatstat.data
## Loading required package: nlme
##
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
##
## collapse
## The following object is masked from 'package:raster':
##
## getData
## Loading required package: rpart
##
## spatstat 1.64-1 (nickname: 'Help you I can, yes!')
## For an introduction to spatstat, type 'beginner'
##
## Note: spatstat version 1.64-1 is out of date by more than 6 months; we recommend upgrading to the latest version.
##
## Attaching package: 'spatstat'
## The following object is masked from 'package:colorspace':
##
## coords
## The following objects are masked from 'package:psych':
##
## reflect, rescale
## The following objects are masked from 'package:raster':
##
## area, rotate, shift
mvest <- st_as_sf(vest, coords = c("Lon", "Lat"), crs = 4326)
str(mvest)
## Classes 'sf' and 'data.frame': 292 obs. of 19 variables:
## $ Uke : int 30 30 30 30 30 30 30 30 30 30 ...
## $ Ã.r : int 2020 2020 2020 2020 2020 2020 2020 2020 2020 2020 ...
## $ Lok_NB : int 30196 15196 10331 12067 12982 10298 11756 12086 13227 10319 ...
## $ Lokalitetsnavn : chr "Ã…dnekvamme" "Aga Ã\230" "Ã…kre" "Aldalen" ...
## $ Lice : num NA 0.36 0.02 0.13 NA 0.18 NA 0.06 NA NA ...
## $ Lus.i.bevegelige.stadier: num NA 7.36 0.55 0.08 NA 0.32 NA 2.37 NA NA ...
## $ Fastsittende.lus : num NA 1.73 0.56 0.18 NA 0.28 NA 0 NA NA ...
## $ Brakklagt : chr "Ja" "Nei" "Nei" "Nei" ...
## $ Telt : chr "Nei" "Ja" "Ja" "Ja" ...
## $ Kommunenummer : int 4634 4613 4617 4624 4645 4615 4632 4617 4624 4602 ...
## $ Kommune : chr "MASFJORDEN" "BÃ\230MLO" "KVINNHERAD" "BJÃ\230RNAFJORDEN" ...
## $ Fylkesnummer : int 46 46 46 46 46 46 46 46 46 46 ...
## $ Fylke : chr "Vestland" "Vestland" "Vestland" "Vestland" ...
## $ Lusegrense.uke : chr "0.5" "0.5" "0.5" "0.5" ...
## $ Over.lusegrense.uke : chr "" "Nei" "Nei" "Nei" ...
## $ Temp : num NA 14.7 16.4 12.4 NA ...
## $ ProduksjonsomrÃ.deId : int 4 3 3 3 4 3 4 3 3 4 ...
## $ ProduksjonsomrÃ.de : chr "Nordhordland til Stadt" "Karmøy til Sotra" "Karmøy til Sotra" "Karmøy til Sotra" ...
## $ geometry :sfc_POINT of length 292; first list element: 'XY' num 5.36 60.83
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "names")= chr [1:18] "Uke" "Ã.r" "Lok_NB" "Lokalitetsnavn" ...
#points$mark <- factor(mvest$Lice)
#points.ppp <- as.ppp(mvest)
#points.ppp$window <- as.owin(poly)
vest2 <- vest[vest$Telt=="JA",]
library(gstat)
##
## Attaching package: 'gstat'
## The following object is masked from 'package:spatstat':
##
## idw
#hscat(Lice~Temp , data=vest2)
g <- gstat(id="Lok_NB", formula=Lice ~ Temp, data=vest2)
str(g)
## List of 3
## $ data :List of 1
## ..$ Lok_NB:List of 15
## .. ..$ formula :Class 'formula' language Lice ~ Temp
## .. .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. ..$ data :'data.frame': 0 obs. of 20 variables:
## .. .. ..$ Uke : int(0)
## .. .. ..$ Ã.r : int(0)
## .. .. ..$ Lok_NB : int(0)
## .. .. ..$ Lokalitetsnavn : chr(0)
## .. .. ..$ Lice : num(0)
## .. .. ..$ Lus.i.bevegelige.stadier: num(0)
## .. .. ..$ Fastsittende.lus : num(0)
## .. .. ..$ Brakklagt : chr(0)
## .. .. ..$ Telt : chr(0)
## .. .. ..$ Kommunenummer : int(0)
## .. .. ..$ Kommune : chr(0)
## .. .. ..$ Fylkesnummer : int(0)
## .. .. ..$ Fylke : chr(0)
## .. .. ..$ Lat : num(0)
## .. .. ..$ Lon : num(0)
## .. .. ..$ Lusegrense.uke : chr(0)
## .. .. ..$ Over.lusegrense.uke : chr(0)
## .. .. ..$ Temp : num(0)
## .. .. ..$ ProduksjonsomrÃ.deId : int(0)
## .. .. ..$ ProduksjonsomrÃ.de : chr(0)
## .. ..$ has.intercept: int 1
## .. ..$ beta : num(0)
## .. ..$ nmax : num Inf
## .. ..$ nmin : num 0
## .. ..$ omax : num 0
## .. ..$ maxdist : num Inf
## .. ..$ force : logi FALSE
## .. ..$ dummy : logi FALSE
## .. ..$ vfn : int 1
## .. ..$ weights : NULL
## .. ..$ degree : num 0
## .. ..$ vdist : logi FALSE
## .. ..$ lambda : num 1
## $ model: list()
## $ call : language gstat(id = "Lok_NB", formula = Lice ~ Temp, data = vest2)
## - attr(*, "class")= chr [1:2] "gstat" "list"
#evgm <- variogram(g)
#evgm <- variogram(g, cutoff=100, width=25)
#revgm <- variogram(g, cutoff=100, width=25, cressie=TRUE)
#cevgm <- variogram(g, cutoff=100, width=25, cloud=TRUE)